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We create a non-Hermitian quantum optimization algorithm to find the ground state of an Ising 
model with up to 1024 spins (qubits). Our approach leads to significant reduction of the annealing 
time. Analytic and numerical results demonstrate that the total annealing time is proportional to 
IniV, where N is the number of spins. This encouraging result is important for the rapid solution 
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I. INTRODUCTION 

Quantum annealing (QA) algorithms can be useful for 
solving many hard problems related to finding the global 
minimum of multi- valued functions, cost functions, opti- 
mal configurations of complex networks, and the ground 
states of the corresponding Hamiltonians. The main idea 
of QA is to utilize the collective quantum tunneling elects 
enabling complex systems to tunnel during the slow time- 
evolution from local minima to the global ground state. 
Although many useful results have been obtained in this 
field, many problems still need to be resolved. The main 
challenge is to accelerate the speed of QA algorithms so 
that the solution time grows, not exponentially, but poly- 
nomially with the size of the problem (l|-[To|. 

There are many approaches to finding the ground state 
of the Hamiltonian, Ho, using QA algorithms (J-Q. One 
of these approaches, is based on introducing a time- 
dependent Hamiltonian, H(t) = Ho + T(t)Hi, m which 
Ho is the Hamiltonian to be optimized, Hi is an aux- 
iliary ("initial") Hamiltonian and [Ho, Hi] ^ 0. The 
T(t)Hi term provides the non-trivial quantum dynamics 
required during annealing. The external time-dependent 
field, T(t), is a control parameter that decreases from a 
large enough value to zero during the evolution. The 
ground state of Hi is the initial state. If T(t) decreases 
sufficiently slowly, the adiabatic theorem guarantees find- 
ing the ground state of the main Hamiltonian, Ho, at the 
end of computation. 

In order to reach the desired ground state, QA algo- 
rithms usually require the existence of a finite gap be- 
tween the ground state and first excited state. However, 



'Electronic address: nesterov@cencar.udg.mx 
t Electronic address: juancarlosbeas@gmail.com 
t Electronic address: gpb@lanl.gov 



in typical cases the minimal gap, g m , is exponentially 
small in the number of qubits, m. For instance, in the 
commonly used quantum optimization m-qubit models 
the minimal energy gap is g m ~ 2~ m / 2 0, 0, ITll4l3| . 
This causes the annealing time to increase exponentially. 
Then, the problem appears on how to speed up the per- 
formance of the QA algorithms and, at the same time, 
find the ground state of Ho ■ 

In [l4| we proposed an alternate adiabatic quantum op- 
timization algorithm based on the non-Hermitian quan- 
tum mechanics. We showed that coupling the system 
to a non-Hermitian Hi produces an effective level repul- 
sion for the total Hamiltonian. This effect leads to an 
adiabatic theory without the usual gap condition that 
rapidly determines the low lying states of Ho, including 
its ground state. 

Recently, we applied this non-Hermitian quantum an- 
nealing (NQA) algorithm to Grover's problem of finding 
a marked item in an unsorted database (l~5| . We showed 
that the energy gap between the ground and the first ex- 
cited state depends on the chosen relaxation parameters, 
and is not exponentially small. The chosen relaxation pa- 
rameters produced a search time that was proportional 
to the logarithm of the number of qubits. 

In this paper we apply NQA to study quench dynamics 
in a lossy 1-dimensional Ising spin chain in a transverse 
magnetic field. We assume that the dissipation vanishes 
at the end of evolution, so that after the annealing is 
finished the system is governed by the Hermitian Hamil- 
tonian. 

We show that dissipation significantly increases the 
probability for the system to remain in the ground state. 
In particular, comparison with the results of the Hermi- 
tian QA reveals that the NQA leads to the ground state 
of Ho with much larger probability, if we use the same 
annealing scheme. We show that NQA has complexity 
of order, hiN, where N is the number of spins. This 
is much better than the quantum Hermitian adiabatic 
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algorithm yielding the complexity of order N 2 . 

This paper is organized as follows. In Sec. II, the 
generic adiabatic quantum optimization (AQO) algo- 
rithm based on non-Hermitian QA is discussed. In Sec. 
Ill, we introduce a lossy 1-dimensional Ising system in 
a transverse magnetic field governed by a non-Hermitian 
Hamiltonian. In Sec. IV, we study the quench dynamics 
of our system using both analytic and numerical meth- 
ods. We conclude in Sec. V with a discussion of our 
results. In the Appendices we present some technical de- 
tails. 



II. NON-HERMITIAN QUANTUM 
ANNEALING: PRELIMINARIES 

The generic adiabatic quantum optimization problem 
based on the QA algorithm can be formulated as follows 
Q . Let T-Lq be the Hamiltonian whose ground state is to 
be found, and %i be an auxiliary "initial" Hamiltonian. 
We consider the following time-dependent Hamiltonian: 



Hr(t) =H Q +g(t)H 1 



(1) 



where [Hq^'Hi] ^ 0. The function, g{t), is monotonic and 
satisfies the relation, g(r) = 0. It is assumed that "Hi is 
dominant at the initial time t = 0, and, since g(r) = 0, 
the H T {t) — r T~io as t — y t. 

The evolution of the system is determined by the 
Schrodinger equation: 



(2) 



The initial conditions are imposed as follows: 1^(0)) = 
\^ g ), where |\&(0)) is the ground state Hi. The adiabatic 
theorem guarantees that the initial state, \^ g ), evolves 
into the final state of |\1/ s (t)), which is the ground state of 
the Hamiltonian, Ho, as long as the instantaneous ground 
state of 'H T (t) does not become degenerate at any time. 
The validity of the adiabatic theorem requires 



E 



{il) m \dU T I 'dt\ip n 



E„ 



< 1. 



(3) 



This restriction is violated near the degeneracies in which 
the eigenvalues coalesce. In the common case of double 
degeneracy with two linearly independent eigenvectors, 
the energy surfaces form the sheets of a double cone. The 
apex of the cones is called the "diabolic point" , and since, 
for a generic Hermitian Hamiltonian, the co-dimension 
of the diabolic point is three, it can be characterized by 
three parameters (IH [I?} ■ For the quantum optimization 
governed by the Hamiltonian of Eq. (fT]) , the requirement 
([3]) can be written as [H H| 



T > T 



max\(i/j e \dH T /ds\ijjg}\ 



mfn I E P 



E g \ 2 



(4) 



where s — t/r, and E e is the energy of the first excited 
state, \4>e)- Thus, if at the time, r c < t, the gap, AE = 



| E e — Eg I , is small enough, the time required to pass from 
the initial state to the final state becomes very large, and 
the AQO loses its advantage over thermal annealing. 

Recently [14j . we proposed a generic non-Hermitian 
adiabatic quantum optimization. Here we consider a par- 
ticular case of the NQA. Let Wo be a Hamiltonian whose 
ground state is to be found, and let, 



Ui(t) = g(l - t/r)Hi - iS(l - t/T)H 2 



(5) 



be the non-Hermitian auxiliary "initial" Hamiltonian. 
Consider the following time-dependent Hamiltonian: 



H T (t)=H +M 1 (t), 



(6) 



where [Ho,Hi(0)] 7^ 0. We impose the initial conditions 
as follows: |*(0)) = $> g ), so that Hi\ip g ) = E g \-ip g ) with 
Eg being the energy of the ground state of the auxiliary 
non-Hermitian Hamiltonian Hi- At the end of evolution 
the total Hamiltonian % t (t) = 'Hqj and the adiabatic 
theorem provides that the final state be the ground state 
of %o i if t ne evolution was slow enough. 

We denote by \^ n {t)) and {tp n (t)\ the right/left 
instantaneous eigenvectors of the total Hamiltonian: 

nr{t)\Mt)) = E n (t)\^ n ), (Mt)\H T (t) = (Mt)\E n (t). 

We assume that these eigenvectors form a bi-orthonormal 
basis, (V> m |^„) = S mn [3. 

For the non-Hermitian quantum optimization problem 
governed by the Hamiltonian the validity of adiabatic 
approximation requires 



T > 



ma,x\(ip e \H T (t)\ip g }\ 



(7) 



mm\E e (t)-Eg(t)\ 2 ' 

where the "dot" denotes the derivative with respect to 
the dimensionless time, s — t/r. This restriction is vio- 
lated near the ground state degeneracy, where complex 
energy levels cross. The point of degeneracy is known 
as the exceptional point, and it is characterized by a co- 
alescence of eigenvalues and their corresponding eigen- 
vectors, as well. Therefore, studying the behavior of the 
system in the vicinity of the exceptional point requires a 
special care [l9l - l2"ll | . 

In the vicinity of the level crossing point, only 
the two-dimensional Jordan block, related to the level 
crossing, makes the most considerable contribution to 
the quantum evolution. Then, the multi-dimentional- 
dimensional problem can be described by the effective 
two-dimensional Hamiltonian, acting in the subspace 
spanned by the ground state and the first excited state 
of the total non-Hermitian Hamiltonian, H T 141. 



III. DESCRIPTION OF THE MODEL 

We consider the 1-dimensional Ising model in a trans- 
verse magnetic field with dissipation governed by the fol- 
lowing non-Hcrmitian Hamiltonian: 



H 



J 



N 

E 

n=l 



9°r, 



'n u n+l 



i25cr„ at) 



(8) 
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with periodic boundary condition, ctn+i — &i- The ex- 
ternal magnetic held is associated with the parameter, g, 
the spontaneous decay is described by the parameter, S, 
and = ( a n ^ ^n)/^ are the spin raising and lowering 
operators. 

We apply the standard Jordan- Wigner transformation, 
and the following procedure outlined in 22-25[, 



— 2r t r 


(9) 




(10) 


2c m c m ), 


(11) 



in which c„ are fermionic operators that satisfy anti- 
commutation relations: {cJ^Cn} = 5 mn and {c m ,c n } = 
{<4,ci} = 0. Then, we obtain H = P+H+P+ + 
P-H-P-, in which 



i N 



(12) 



denote the projectors onto the subspaces with even " + ' 
and odd " — " numbers of quasiparticles [H, HH , and 



H ± = -— ^2 {cjc n+1 + c n+1 c n + g + iS 



N 



^9 c n c n + C n+lCn + c n C n+l) > (13) 



in which g = g + i8. 

The Hamiltonian, H~ , is related to c„'s with the pe- 
riodic boundary conditions, cn+i — c%, while in H + the 
operators, c„, obey the following ('antiperiodic') bound- 
ary conditions: cyv+i = —c\. Since the parity of quasi- 
particles is conserved, one can consider only either H + 
or H~ . Further we consider only quasiparticles with the 
even parity. 

Applying the Fourier transformations with the an- 
tiperiodic boundary condition, cn+i = — Ci, 



E 



c k e 



Akn 



(14) 



in which k — ±ir/N, . . . , ±(N — l)ir/N is half-integer 
quasimomentum, we obtain, H + = ^ fe Hk, where 



Hk = |( 2 (3- cos k ) c l°k ~ 9 



+ sin k(c\c\ k + c-kCk)). (15) 
The Hamiltonian, H + , can be diagonalized by using 



the generalized Bogoliubov transformation: 



"k " ~k f 

Cfe = cos — bk + sin -y- V_ 



k ' 



where 



4 = cos y 6 fc + Sm — h -k' 

, #fc . . Ok + 

Ofe = cos yc fc + sin — c^ fc , 

K = cos y c fe + sin y c -fe> 



COS fc 



sin y fc 



y/g 2 — 2<?cos fc+T 
sin fc 

\/<7 2 — 2g cos fc + 1 



(16) 
(17) 
(18) 
(19) 

(20) 
(21) 



with 9k being a complex angle. 

There are two eigenstates for each fc, 



M*)> = (£4 



(u + (fc)| = (cos^,shA, (22) 



-Ml = (- sm y, cos y), 



(23) 



with the complex energies, s±(k) = —So±ek, where £0 = 
Jcosfc + iJS, and £k = J^/ g 2 — 2g cos k + 1. Here we 
denote by |u±(fc)) ((u±(fc)|) the right (left) eigenvectors. 

With help of Eqs. (fT6|) - ([19]) we obtain the diago- 
nalized Hamiltonian as a sum of quasiparticles with half- 
integer quasimomenta, 



ft fc 

E^+E^( 6 ^- 6 -^-fe)- (24) 



fc>0 



fc>0 



Its spectrum contains only states with even number of 
quasiparticles. 

The ground state of the Hamiltonian (|2"4")l is a state \tp g ) 
annihilated by all quasiparticles annihilation operators, 
bk, so that, bk\ip g ) = 0. One can show that the ground 
state can be written as a product of qubit-like states: 



l^>=0 fe (cos^|O) fe |O)_ 



fc-sinf |l) fc |l) 



,(25) 



<V> 9 1 = ® k [ cos % <0| fc <0| _ fc - sin ^ <l| fc <1 1 - fe J , (26) 

where, |0)fc, is the vacuum state of the mode, Cfc, and |l)fc 
is the excited state: |l)fc = cj.|0)fc. 

Since for each fc, the ground state lies into the two- 
dimensional Hilbert space spanned by |0)fc|0)_fc and 
|l)fc|l)_fc, it is sufficient to project Hk on this subspace. 
For a given value of fc, both of these states can be repre- 
sented as a point on the complex two-dimensional sphere, 
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S 2 , with coordinates, (6k, ip). In this subspace the Hamil- 
tonian, Hk, takes the form 

For \g\ 3> 1, the ground state is paramagnetic with all 
spins oriented along the x axis, and from Eq. (|20[) we ob- 
tain cos 9k — y 1 as |<7| — y oo. Thus, the south pole of the 
complex Bloch sphere corresponds to the paramagnetic 
ground state. On the other hand, when \g\ <C 1 there are 
two degenerate ferromagnetic ground states with all spins 
polarized either up or down along the z axis. The real 
part of the complex energy reaches its minimum at the 
point defined by cos 9k = — 1, and, hence, the north pole 
of the complex sphere is related to the pure ferromag- 
netic ground state with the broken symmetry in which 
all spins have orientation either up or down. However, 
in the thermodynamic limit the system passing through 
the critical point ends in a superposition of up and down 
states with finite domains of spins separated by kinks 

The ground state energy is given by, 

E gs =-Ne -J2 £ k> ( 28 ) 

fc>0 

and in the thermodynamic limit (N — ¥ oo) the ground 
state energy per spin can be written as 



e gs =-£ -- / y/g 2 -2<? cos (29) 

T JO 

where the angle § = k as N — > oo. This integral can 
be written in terms of a complete elliptic integral of the 
second kind, 

£gs = -US - 2J( ^ + 1} E(2Vl/(g + 1)). (30) 

Fig. [T] shows in the thermodynamic limit the absolute 
value of the difference between the two eigenvalues of the 
Hamiltonian (|38[) . 

|Ae| = 2J|V<? 2 -2<7costf + l|, (31) 

as functions of g and $. As one can see, the energy gap 
vanishes at the critical point, 



d c = arccos y/l - S 2 , (32) 
g c = y/1 - S 2 . (33) 

The difference between the Hermitian QA and non- 
Hermitian QA is that, while in the first case the gap van- 
ishes for long wavelength modes (i? c = 0), in the second 
case the minimal gap shifts to short wavelength modes 
[§ c = arccos \/l — 5 2 ). 




FIG. 1: (Color online) Absolute value of difference, Ae, of the 
eigenvalues of the Hamiltonian (|27|l as functions of g and ■& 
in the thermodynamic limit. Left panel: S — 0. Right panel: 
S = 0.25. 

IV. QUENCH DYNAMICS 

In this section, we consider the NQA for the time- 
dependent Hamiltonian of Eq.® written as 

Hr(t)=U +Hi(t), (34) 

where 

J N 

«o = ~£«+i' (35) 

71=1 
J N 

= - 2 E + im)°n<)- (36) 

n=l 

We start with the ground state of the auxiliar Hamil- 
tonian, %i(0), as the initial state, which is a "paramag- 
netic" with all spins oriented along the x axis. For g ^> 1 
the Hamiltonian, H T (0), is dominated by Hi(0), and the 
ground state of the total Hamiltonian, "H T , is determined 
by the ground state of 'Hi(O). The Hi term causes quan- 
tum tunneling between the eigenstates of the Hamilto- 
nian, Ho. At the end of NQA we obtain, H t (t) = Ha- 
lt the quench is slowly enough, the adiabatic theorem 
guarantees reaching the ground state of the main Hamil- 
tonian, Hq, at the end of computation. 

As shown in the Sec. Ill, the total Hamiltonian, 
H T (t), in the momentum representation splits into a sum 
of independent terms, H T {t) = 5Z fe Each Hk 
acts in the two-dimensional Hilbert space spanned by 
\ki) = |l) fc |l)-fc and \k ) = |0)fc|0)_ fe . The wave function 
can be written as, \ip) — Y\ k \ipk), where 

\M*))=co(k,t)\k ) + c 1 (k,t)\k 1 ). (37) 

Choosing the basis as, |fci) = (^j and |fco) = (^j , 

we find that the Hamiltonian, Hk{t), projected on this 
two-dimensional subspace takes the form 

wo-^oi+^^r* -«*«*) - 

(38) 



5 



where £o(t) = J cosk + iJS(t) andg(i) = g(t)+iS(t) [48 
Further we assume linear dependence of g(t) on time: 



m = 



l(r-t), 0<t<r 
0, t > T 



(39) 



where 7 = (g + iS)/r, and g, (5 are real parameters. 



A. Diabatic basis 

The general wave functions, \ip k ) and (ip k \, satisfy the 
Schrodinger equation and its adjoint equation 



W*(i)|Vfc), (40) 

: mU k (t). (41) 

Presenting ^(i)) as a linear superposition, 

|lk(*)> = (ufc(*)|*0> + «fc(*)|fci»e^ eoW *, (42) 
and inserting expression (|4"2"j) into Eq. (|40p . we obtain 

mfc = j( — (5 — cos fc) Ufe + sin kvk), (43) 



W fe = J( sin k u k + (g - cos fc) v k ) . 



(44) 



The solution can be written in terms of the parabolic 
cylinder functions, D^i Vk (±z). (For details see Appendix 
A): 

U k (z k ) = B k D- iUh (z k ) - iVu^A k D it/h -i(iZk), (45) 
Vfe(z fe ) = A k D iVk (iz k ) - \fw k ~B k D- lUk ^ 1 (z k ), (46) 

in which we introduced new functions: u k {t) — U(z k ), 
v k (t) = V(z k ) and 



z k (t) 



t/4 



cos k 



v k 



J sin k 

27 



(47) 
(48) 



The constants, and Bfc in Eqs. (|45l) . (T4"6"l) . are deter- 
mined by the initial conditions. 

In what follows we assume that the evolution of the 
spin chain starts at t = in the 'ground' state. This im- 
plies that for each k the evolution of the corresponding 
two-level system (TLS) starts from the state, |?/>(0)) = 
|fco). Then, the following initial conditions should be im- 
posed: Ufe(0) = 1 and v k (0) = 0. The related bound- 
ary condition are, z k (0) = e"^yj2 Jpf^yr — cosfc). 
Using these conditions, we obtain the solution of the 
Schrodinger equation as follows. (See Appendix A for 
details.): 



U k (z k ) 
V k {z k ) 



B k D^ iVh (z k ), 



nj7^.B k D- il/h -i(z k )), 



(49) 
(50) 



where B k = e™«' 2 D tVk (iz k (0)). 

It is assumed that a quantum measurement will deter- 
mine the state of the quantum system at t > r, when the 
external field, g(t) = 0. We denote the final state of the 
system, at t = r, as \ip T ) — J\ k \ip k (r)). The probability, 
P n (k), of finding the TLS in a given state, \k n ) (n = 0, 1), 
can be written as 



Pn(k) 



l(fcn|^)| 2 



(51) 



Since for non-Hermitian systems the norm of the wave- 
function is not conserved, we define the (intrinsic) prob- 
ability of transition, Ifco) — > |fcx), as 



Pk(t) 



\Mt)\ 2 



(52) 



Using the functions, U(z k ) and V(z k ), we recast (|52|) as 

1 



P k {t) = 



1 + 



\D. 



.(^(*))l s 



(53) 



'Tv k ~D- iVk -i{z k (t))\ 2 



To calculate P k (t) at the end of evolution it = r) we 
use asymptotic formulas of the Weber functions. For 
large values of the argument, |zfc(r)| = lv2J/7Cosfc| 
1, one can apply the asymptotic formulas for parabolic 
cylinder functions to estimate the probability of transi- 
tion. For r S> 1 the modulus of this argument is large 
for most fc, except near ka = ±7r/2. 

For wavelength modes with ka < 7r/4, using the 
asymptotic formulas for the Weber functions, we obtain 



D- 



,(%(r)) 



e -^/2 e -^W/2 r(1+in ) 



Inserting (|54j) into Eq. (|53|) . we obtain 

1 



Pk{r) = 



1 



\T{l + iv k ) 
2ir\v k \ 



(54) 



(55) 



3 -7rSRy fc -5Rz^(r) 



For 5 -C g we can approximate T(l + iv) »r(l| ffilv) 
Next, using the relation (27| . 



for real y, we obtain 
Pk(r) = 



?/ sinh iry 



I _ e -1ixlSiv k _|_ e -2^5Ri/ fc -SRz| f (t) ' 



(56) 



(57) 



In the case of the Hermitian QA (<5 = 0),one has 
^tz\(r) = 0. and Eq. ([55)1 leads to the Landau-Zener 
formula [HlH, 



Pk(r) 



1 



-(ttJt/ g) sin 2 



(58) 
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FIG. 2: The probability of transition (remain in the ground 
state), |fco) —> \ki), as a function of the scaled time, s — t/r 
(J = 0.5, g = 10, r = 10 3 , N = 128). Left panel: Hermitian 
QA (5 = 0). Right panel: NQA (S = 0.5). Blue line (k = 
7r/128), red line (k — tt/32), green line (k = tt/8), black line 
(k — 7r/4), orange line (k — tt/2). 



To calculate the transition probability for k ps tt/2 we 
use the expansion for small value of the argument of We- 
ber function, \z k \ <C 1 (see Eq. (|A44|> ). The computation 
yields 



P*/2(t) = 



1 



|r(W 



(59) 



|r(i^)|s 



For the Hermitian QA this gives 



P 



V/2W = - 



tanb.7rz//2 
+ tanh.7rz//2 ' 



(60) 



in which v = Jr/(2g). From Eq. (jST))) it follows that for 
t > g/J the probability P 7r / 2 ('r) ^ 1/2. 

In Fig. [5] we present the results of our numerical sim- 
ulation for N = 128 qubits. As one can see, for long 
wavelength modes with k <C 7r/4 (blue and red lines), 
the NQA shows better performance than the Hermi- 
tian QA. For k = tt/2 the probability of transition is 
Ptt/2( t ) — 1/2 for both schedules, either Hermitian QA 
or non- Hermitian QA (orange line). 



B. Adiabatic basis 

The widely used opinion is that for slowly enough evo- 
lution the only long wavelength modes become excited 
(see e.g. [23j]). However, as it was shown in the Subsec- 
tion A, even for the Hermitian QA, the transition proba- 
bility, Pfc fa 1/2, for k ps 7r/2. Thus, it is not clear what is 
the contribution of the short wavelength modes (k ~ 7r/2) 
to the probability of the whole system to stay in the 
ground state. To respond to this question we consider 
the expansion of the wavefunction, \ip) = J| fc \ipk(t)}, in 
the adiabatic basis, formed by the instantaneous eigen- 
vectors. 



In the adiabatic basis the wavefunction, (^(f)), can 
be written as follows: 



\i/>k(t)) = a k (t)\u-{k,t))+l3k(t)\u + (k,t)). 



(61) 



We assume that the evolution begins from the ground 
state that implies, aifc(O) = 1 and /?fc(0) = 0. At the end 
of evolution at t = t, when g = 0, we have 



|Vfe(r)) = a k (T)\u-(k,T)) + p k (T)\u + (k,r)), 
where 



\' u +( k , T )) = 
\ u -( k , T )) = 



sin ■ 



By presenting, 

a k (t) = a k (t)e- i ti s -^ dt , 
p k {t) = b k {t)e- i K e +^ d \ 



(62) 

(63) 
(64) 



(65) 
(66) 



one can show that the coefficients, a k (t) and b k (t), satisfy 
the following asymptotic conditions: 



o fc (r)=l + 0(l/T), (67) 
b k (t) = o( cxp (2r3^ " e k (z)dz) \ (68) 



where the critical point, z c , lies on the first Stokes line in 
the lower complex line defined as 



J / e k {z)dz < 0. 
'o 



(69) 



Here the critical point, z c , is determined as a solution of 
the equation, e k (z c ) — 0, in the c omp lex plane obtained 
by analytical continuation, t — > z |3CH38| . 

Similarly, if initially the system was in the excited 
state: \ip k (0)) = \u + (k,0)), so that a k (0) = and 
/3fc(0) = 1, the result of integration yields 

b k (r) = l + 0(l/r), (70) 

a k (r) =0^exp(2r3^ ° e k (z)dz^ . (71) 

The intrinsic probability to remain in the ground state 
at the end of the adiabatic evolution is given by 



P a k \r) = 



Mr)| 2 + |/3 fc (r)r 
With help of Eqs. (|65|). (l66t we obtain 

P^ir) = ' 



1 



\ b k(T)\ ACj fr eh(t)dt 



(72) 



(73) 



7 



Prom here it follows that for any, < t < r, the adia- 
batic evolution should be performed along the path cor- 
responding to, SsE k (t) < . In the exact solution, given 
by 



Qk 



(r) = - B k e- 5r ' 2 (D_ iVk (z fc (r)) sh 



+ ^/IU k ~D- Wk - 1 (z k (T))) cos 
3 fc (r) =B k e- 5r / 2 ( D_ lVk (z k (r)) cos 



- s/iUkD - il/k -i(zk(T))) sin - ), 



(74) 



(75) 



it manifests itself in the choice of phase in the argument 
of the Weber function, when we apply the asymptotic 
expansion. 

The probability for the whole system to stay in the 
ground state at the end of the evolution is given by, 
Pgs = n k P k S { T )- F° r l° n S wavelength modes with 
k <C 7r/4, using the asymptotic formulas for the Weber 
functions with the large value of its argument, we find 
that P k s (r) = P k {r), where P k (r) is the transition prob- 

For the Hcrmitian 



ability of spin flip given by Eq. (E 
QA this yields the LZ result 



Pk(r) = 1 



-(it Jt / g) sin 2 k 



l 



-TV 3 jTk 2 /(gN 2 ) 



And in the case of the NQA (for S <C g) we obtain 



Pk{r) 



where 



1 _ g-27rSfti/ fc i g-27rSRi/ fc -fRz£(-r) ' 

3?i/ fc ps (Jr/2 5 )sin 2 fc, 
^zI{t) ps {26 Jt / g 2 ) cos 2 k. 



(76) 



(77) 



(78) 
(79) 



For short wavelength modes, approximately with 
7r/4 < k < 7r/2, employing the large-order asymptotics 
for Weber functions, we obtain Pf s (r) = 1 + 0(1/ y/\v k \). 
(See Appendix A.) 

Our theoretical predictions are confirmed by numerical 
calculations performed for N = 512 qubits. (See Figs. |31 
|4j) One can observe that while short wavelength exci- 
tations are essential at the critical point, at the end of 
evolution their contribution to the transition probability 
from the ground state to the first excited state is negli- 
gible. 



V. PERFORMANCE CHARACTERIZATION OF 
THE QUANTUM ANNEALING 

During the QA the system does not stay always at the 
ground state at all times. At the critical point, the sys- 
tem becomes excited, and its final state is determined by 
the number of defects (kinks). To evaluate the efficiency 
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FIG. 3: (Color online) The probability, P| s , of TLS to re- 
main in the ground state as a function of the scaled time, 
s = t/r, for the Hermitian QA (S = 0, J = 0.5, g = 
10, r = 10 3 , N = 512). From bottom to top: k = 
tt/512, tt/256, tt/128, tt/32, it/16. 
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FIG. 4: (Color online) The probability, P| a , of TLS to remain 
in the ground state as a function of the scaled time, s = t/r, 
for the NQA (5 = 0.5, J = 0.5, g = 10, r = 10 3 , N = 512). 
Blue line (k — 7r/512), red line (k = 7r/256), green line (k = 
7r/128), orange line (k — tt/32). 



of QA one can calculate the number of defects. Then, 
computational time is the time required to achieve the 
number of defects below some acceptable value. 

Following [23|, we define the operator of the number 
of kinks as 

1 N 



fe>0 



The number of kinks is equal to the number of quasipar- 
ticles excited at g = (final state). It is given by 



Af = (Vvl-A^lVO = jYl ((Mr)\n k (r)\^ T ) - E gs ). 



fc>0 



(81) 



This can be recast as follows: JV = E res /J, where E res 
is the residual energy defined as the difference between 
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the solution obtained by Q A and the exact one [!, H, H^] , 

E res = ((Mr)\n k (T)\4, T ) - E gs ) = JN. (82) 

fe>0 

Thus, an alternative (equivalent) way to evaluate the ef- 
ficiency of QA is to calculate E res . Note, that for more 
complicated systems that include disorder, the residual 
energy may not have a simple link with the density of 
defects Hf. 

Using Eq. (j6"2")l , we can calculate the number of kinks 

as 



AA=£(1-Pf(r)), 



(83) 



k>0 



where P{! s {t) is given by Eq. (|72|) . 

In the thermodynamic limit the sum in Eq. (|83p can 
be replaced by integral, and we obtain for the density of 
kinks the following expression: 

n= lim ^ = - f dk{l-P a k s {r)). (84) 

N^oo TV I? Jo 

As was shown in the previous section, during slow evolu- 
tion only long wavelength modes can be excited. So one 
can use the Gaussian distribution by replacing sin k sa k 
and cos k ~ 1 . In the limit a/2 Jr/g ^> 1 , we can employ 
Eqs. (|77|) - ([79)) to calculate the number of kinks as 



1 

n = — 

* Jo 1 



e -2irSR;/ fc _|_ e -2^u k -Mz 2 k (r) ' 



(85) 



Performing the integration with = Jrk 2 /(2g) and 
3^(t) = 25Jr/g 2 , we obtain g| 



-2SrJ/g 2 



" 2 



,1 



where 



710 = 9^ 
Z7T 



(86) 



(87) 



denotes the density of kinks for the Hcrmitian LZ prob- 
lem [23j], and $(x,a, c) is the Lerch transcendent [4l| . 

In Figs. Oiniwe present the results of numerical simu- 
lation for the density of defects. In Fig. [5] the density of 
kinks as a function of 6 is depicted. In Fig. [5] we show 
the dependence of the density of defects as a function 
of the decay parameter, S, and annealing time, r, in the 
thermodynamic limit. This is consistent with the results 
of numerical simulation presented in Ref. [42[ for the 
Hermitian LZ problem. 

The final state of the system is a ferromagnetic state 
with the finite domains of spins (pointed up or down), 
separated by kinks. The magnetization, TVm 2 , is defined 
from the expression for the total energy of the ground 
state: E gs = —NJm 2 z . In the thermodynamic limit, we 
obtain 



m z = lim 

IV->oc 



E 



ys 



JN 



(88) 
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FIG. 5: (Color online) Density of kinks as function of the 
dissipation parameter, 5 (r = 10 3 , J = 0.5, g = 10, TV = 
1024). Blue line: exact result. Red dashed line: asymptotic 
formula of Eq. ((86|) . 
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FIG. 6: (Color online) Density of kinks as function of the 
dissipation parameter S and annealing time r. 



In Fig. [8] the density of magnetization as function of 
the annealing time, r, is depicted. As one can see, even 
moderate dissipation essentially decreases the number of 
defects in the system. 

Due to the symmetry of the problem with respect to 
k — > — k, and since for each k there is independent evo- 
lution, the probability of the whole system to remain in 
the ground state at the end of evolution is the product 



(89) 



fe>0 



In the long wavelength approximation one can take 
into account only k — 7r/TV, and estimate P gs as 



P, 



1 



(90) 



where 3?z 2 (r) = 26 Jr/g 2 and div = t/tq. Here we denote 
t = 2gN 2 /(7r 2 J). 

In Figs. H] and [TU1 the results of numerical simulation 
are demonstrated. As one can see, for the probability, 
P gs « 1, the asymptotic formula of Eq. (|9T)]) is in a 
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FIG. 7: (Color online) Density of kinks as function of the 
dissipation parameter, 5 (J = 0.5, g = 10, r = 10 3 , TV = 
1024). Blue solid shows the exact result. Dashed color lines 
present contribution of the first n-modes. From bottom to 
top: n= 1,16,32,64 . 
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FIG. 8: (Color online) Density of magnetization, m z , as 
function of the annealing time, r. From bottom to top: 
5 = 0, 0.25, 0.5, 1 (J = 0.5, g = 10). 



good agreement with the exact formula (|89j) . We also 
performed numerical simulations to demonstrate that for 
any TV the contribution of the first TV/16 modes yields 
essentially the same result as the exact formula (|89p . (See 
Fig. |9l) We find that even moderate dissipation boosts 
the transition probability. 

For the Hermitian QA (6 = 0), Eq. flSD) yields the 
Landau-Zener formula 12 a, l29l 
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FIG. 9: (Color online) Probability to stay in the ground 
state, P gs , as function of the dissipation parameter, 8 (J = 
0.5, g = 10, t = 10 3 ). From top to bottom: TV = 
64, 128, 256, 512, 1024. Left panel: Solid curves present the 
exact result. Dashed lines present the contribution of the 
first mode (k — n/N). Right panel: Solid curves are the con- 
tribution of all modes. Dashed lines are the contribution of 
the first TV/16 modes. 
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FIG. 10: (Color online) The transition probability, P T , as 
function of a scaled decay rate, 5* = Sro/g 2 , and scaled an- 
nealing time, t* = t/to- 



P gs = l-e- 2 * T ' T °. 



(91) 



From here it follows that P gs rj 1, if r > tq. Thus, the 
computational time for the Hermitian QA should be of 
order TV 2 . 

For the NQA, assuming r <C tq, we obtain 



P, 



1 



1 + ^L e Sr/ 9 > 
2lTT 



From here, in the limit of S — > 0, we obtain 

1 



1 



2ttt 



< 1. 



(92) 



(93) 



The obtained result is expected, as in this case, the time 
of the Hermitian annealing, r, is small with respect to 
the characteristic time, tq: t <C tq. 



Next, assuming 



^-ln-^»l. 
g A Zttt 



we obtain 



Pa, 



2-ITT 



(94) 



(95) 



As one can see, P gs ~ 1, if the conditions of Eq. 



(|94|) are satisfied. From (|94|) we obtain the follow- 
ing rough estimate of the computational time for NQA: 
t fa ( 5 7<5)lnTV. 

To find the annealing time from the exact formula (|90)) 
we impose the following condition on the probability: 
P T = 0.999. Then, we solved numerically Eq. 0- The 
results of numerical calculations are presented in Fig. 1111 
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FIG. 11: (Color online) Annealing time, r, as a function of 
N for NQA (J = 0.5, g = 10). Solid lines present the exact 
result. Dashed lines are the asymptotic formula of Eq. (|96l) . 
From top to bottom: 8 = 0.05, 0.1, 0.25, 0.5, 1. 
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VI. CONCLUSIONS 

Recently, many modifications of quantum annealing al- 
gorithms have been proposed [IHH, H4IS El S3- The 
main objective of these publications is to significantly de- 
crease the time of annealing so that the solutions of hard 
optimization problems could be obtained either by (i) 
combining classical computers with quantum algorithms 
or by (ii) building real quantum computers. One of the 
very popular test models is the Ising spin system that is 
also useful for practical purposes. In this case, the quan- 
tum annealing algorithms are used to find the ground 
state of this system. The approach presented in this pa- 
per, is contained in item (i). We choose an auxiliary 
Hamiltonian in such a way that the total Hamiltonian is 
non-Hermitian. This allows us to shift the minimal gap 
in the energy spectrum in the complex plane, and signifi- 
cantly reduce the time required to find the ground state. 
Our approach leads to the annealing time ~ InN, where 
N is the number of spins, which is much less than the 
time of Hermitian annealing (~ N 2 ) for the same prob- 
lem. But many serious problems still remain to be consid- 
ered. One of them is the application of this dissipative ap- 
proach to more complicated Ising models with frustrated 
interactions, when the ground state can be strongly in- 
homogeneous. This research is now in progress. 



FIG. 12: (Color online) Annealing time, r, as function of 
N. Left panel. Hermitian QA, 5 = (J = 0.5, g = 10). 
Right panel. Non-Hermitian QA, from top to bottom: 8 = 
0.05,0.1,0.25,0.5,1. 



We find that the best fit of the asymptotic expression to 
the exact result is given by the following asymptotic for- 
mula (see Fig. ITTj) : 
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~ ^ i N rtwrt 
T ~ ~~S~ tt' Appendix A: Exact solution of the Non-Hermitian 

Landau-Zener problem 



In Fig. [T2j the annealing time, t, as a function of 
the number of spins is depicted for NQA (left panel) and 
Hermitian QA (right panel). The comparison shows that 
for large number of spins (N ~ 1000) and for 6 > 0.25, 
the annealing time of NQA is ps 10 3 times smaller than 
for Hermitian QA. 

The obtained results indicate that the characteristic 
time of non-Hermitian annealing, even for small but finite 
6 0, is defined not only by the number of spins, N (as 
in Hermitian annealing), but mainly by the dissipation 
rate, 5. (See Figs. El [TO] and Eq. (p]t.) Thus, the non- 
Hermitian quantum annealing has complexity of order 
In AT, which is much better than the quantum Hermitian 
(global) adiabatic algorithm. Also, this complexity is 
certainly better than one of the adiabatic local annealing 
algorithm which has a total running time of order N (43j . 



The non-Hermitian Hamiltonian, Hk(t), projected on 
the two-dimensional subspace spanned by \k%) = ( J J 

and |fco) = f i ) i takes the form 

(Al) 

where £o(i) = Jcosfc + iJ5{t) and g(t) = g(t) + iS(t). 
We assume a linear dependence of the function, g(t), on 
time: 

where, 7 = (g + iS)/r, and g, S are real parameters. 



11 



The general wave functions, \ip k ) and (ip k \, satisfy the 
Schrodinger equation and its adjoint equation 



i-\rPk) = H k (t)\^ k ), 



Presenting \tp k {t)) as a linear superposition 

\Mt)) = (u k (t)\k ) + Vk {t)\k l )ys^ d \ 

and inserting (|A5|) into Eq. (|A3|) . we obtain 

iu k = J{ — (g — cos k) u k + sin kvk), 
iv k = J{ sin kiik + (g — cos k)v k ). 



(A3) 
(A4) 

(A5) 



(A6) 
(A7) 



the parabolic cylinder functions with — 37r/4 < argz^ < 
37r/4, we obtain 

A k = + O(l/|^(0)|), (A16) 

B k = (z k {Q)r k e zl(0)/i - (A17) 

Similar consideration of the adjoint Schrodinger equa- 
tion with the wavefunction 



{i>k\ = (u k (t)(ko\+v k (t)(ki\)e-^ £o{t)dt , 



yields 



iu k = — J( — (g — cos k) u k + sin k v k ) 
iv k = — j( sin ku k + (g — cos k)v k ). 



(A18) 

(A19) 
(A20) 



For the functions, u k {t) = U k (z k ) and v k (t) = V k (z k ), we 
obtain 



Let z fe (t) = e 47r / 4 ^2J/7( 7 (r - t) - cosk) be a new 
variable. Then, for new functions, u k (t) = U k (z k ) and 
v k (t) — Vk(zk), we write Eqs. (|A6|) . (IA7|) in the standard 
Landau-Zener form [28l Sal 



-j—Uk = ^-Uk - V^kV k , 
dz k 2 

■3— Vfc = — — V k - \fiv k U k . 
dz k 2 



dz. 



-Uk 



Zk 



Vk 



d z k 

V k = —V k 

dz k 2 



IWkVk-, 



(A8) 
(A9) 



From here it follows 
dzf 



"k 

dz 



(5 



2 V k 



+ ^+iv k )U k =0, 



iv k )V k = 0. 



where v k = Jsin 2 fc/2 7 , and the complex 'time' z k The solutions are given by 
runs from z k (0) = e l7T / 4 y / 2 JpyijT — cos A:) to z k (r) = 



-e 47r / 4 v/2j77Cosfc. 

From Eqs. (|A8J) . (|A9|) we obtain the second order We- 
ber equation 



U k (z k ) 

V k (z k ) : 



A k D- iVk ^i(z k ) + B k D iVk (iz k ), 
1 



(A21) 
(A22) 

(A23) 
(A24) 

(A25) 



>iv k 



--A k D_ iUk (z k ) + is/E> k B k D iUk -i(iz k ), 

(A26) 



^ fc )[/ fc = 0, (A10) 
(All) 



where 



Solution of the Weber's equation is given by the parabolic 
cylinder functions _D_.;„ fc (±z), Di„ k -i(±iz). 

We obtain the solutions of Eqs. (|A81 IA9j) in the form 



(A27) 
(A28) 

(A29) 
(A30) 



U k {z k ) = B k D_ iUk {z k ) - iy/iu k ~A k D iVh -i(iz k ), (A12) 
V k (z k ) = A k D iVk (iz k ) - y/iv k ~B k D- iVk -i(z k )), (A13) 

where the constants, A k and B k , are determined from 
the initial conditions. 

If the evolution of TLS starts at to = in the 'ground' 
state, \ip(0)) = \k ), the following initial conditions 
should be imposed: u k (0) — 1 and v k (0) = 0. Using 
the identity $EMj , 

we obtain 

A k = ^e^^D^^z^O)), (A14) 
B k = e™»/ 2 D Wk (iz k {Q)). (A15) 

We assume further that r g/J. This implies 
|^fc(0)| 3> 1. Then, applying the asymptotic formulas for 



A k = ^e™ fc / 2 A„ fc -i(^fc(0)), 

B k = e™»/ 2 D^ lUk (z k (0)). 

For |zfe(0)| ^> 1 we obtain 

A k =Q + O(l/|%(0)|), 

B k = e™"/ 2 {z k {Q))- iVk e" 2 ^°)/ 4 . 

Finally, the straightforward computation shows that the 
obtained solutions satisfy the normalization condition 

(^fc(t)lV'fc(*)) = i- 

The solutions of the Schrodinger equations for 
I (0)| >• 1 and the initial conditions: u k (0) = u k (0) = 1, 
«fc(0) = v k (0) = 0, are given by 

U k (z k ) = B k D_ iVk (z k ), (A31) 

V k (z k ) = -y/ivkB k D- iUh - 1 (z k ), (A32) 

U k (z k ) = B k D iVk (iz k ), (A33) 

V k (z k ) = i\pw~ k B k D iVk _ x (iz k ), (A34) 



where 



B k = (z k (0)r* e^'\ 



B k 



(*fc(0)) 



-™k _-2;(0)/4 



(A35) 
(A36) 
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1. Some important properties of the Weber 
functions 



The parabolic cylinder functions, D_j„ fc (±2), 
Di Uk -i{±iz), being solution of the linear differen- 
tial equation 



In particular, for — hir/A < &rgz k < — 7r/4, this yields 
D^ k (z k )=z-^e-^\l + 0(\zl\- 1 )) 



iV2vr 
T(iv k ) 



e-* v *z% h - 1 <?>''\\ + 0\zl\- 1 ) 



(A47) 



d 2 U (\ 



- + U- — }U = 



dz 2 \2 4 



(A37) 



b. Large-order asymptotics 



satisfy the following derivative and recurrence relations 

- d ^{ eZl/4D -^)) =iVke z * /4 D^ k ^(z k ), (A38) 
-^- k [ e ~ zl/4D -^( z ^) =-e- 4/4 D- iUk+1 (z k ). (A39) 

D-i Uk+ i(zk) - z k D_ iVk (zk) + v k D- iVk -i(z k ) = 0. 

(A40) 

The Wronskian for these solutions is given by 
D- t v k (z k )-^-D- lUk (~z k ) 

D_ iUk (-z k )-^-D_ iUk (z k ) = (A41) 
dz k T(iv k ) 

D- iVh {z k )—D iVk+ x{izk) 

- D iVk+1 (iz k )J-D- iVk (z k ) = -is™*' 2 . (A42) 
dz k 



For large-order value of the Weber functions with a 
phase of argument | argz^l < ir/2 the leading terms are 

sua 



D. ivk (z k ) ~ cos % e «*/4-*> (l + ) . (A48) 

D_ iVk _x{z k ) ' ' ' 



^ e W4-wy( 1 + (D 



where 
and 



2 

1 . 
— — sin 

'iu k 2 



i7r/4 



(A49) 



cos6*fc = 



z fc e 



t/4 



(A50) 



(A51) 



Using Eqs. (jA38|) - (|A42)l . we obtain 

D- iUk {z k )D iUk (iz k ) + is k D_ i „ k _ 1 (z k )D i „ k _ 1 (iz k ) = e~™ k/2 

(A43) 

For small value of the argument one can use the power- 
series expansion of Weber function yielding 

D_ Wk (z k )=2-^/ 2 ^ as|z fc |->0. (A44) 

r (2 + -±) 

a. Asymptotic expansion for large value of argument 

For large value of argument, \z k \ ^> 1, and for 
| arg z k \ < 37r/4 the following asymptotic expansion is 
valid [|l| 

D_ lUk {z k ) = z-^e- z *l\l + 0{\z 2 \- 1 )) (A45) 

To find the asymptotics of the Weber functions for other 
values of its argument one can use the relations: 



Appendix B: Equation of motion 

We consider a two-level system governed by the non- 
Hermitian Hamiltonian, 'Heff, written as 



*U = Yi + in(t)-«7, 



(Bl) 



where Hit) — ((l x (t),D, y (t),Cl z (t)) is the complex vector 
and A = A - iT, where T = (r + r x )/2. The qubit 
states \u(t)) and (u(t)\ satisfy the Schrodinger equation: 



=n eS Ht)), 



(B2) 
(B3) 



D- iVk {z k ) = e-™<°D- iVk {-z k ) 



iV2ir 
T(w k ) 



e-^l 2 B iVk _ x {iz k ). 
(A46) 



Employing Eqs. ([B2|) . (|B3|) . we find that the Bloch vec- 
tor, n(t) = (u(t)\cr\u(t)) , satisfies the following general- 
ized Bloch equation (GBE): 

— =-Tn + n^fl(t) + 3tn(t) x n, (B4) 
dt 



^ = -r» + 9M(t)-n, 
dt 



(B5) 
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where n = (u(t)\u(t)) = J n 2 x + n 2 + n 2 z . PooW- Substituting these expressions for n z (t) and n(f) 

V - into Eqs. |B6} - (IB"9l) . we obtain 
Denoting the real part of the complex vector $7 as 

!Jff2 = (p. x ,Q, y ,n z ) and its imaginary part as 3f2 = 

(A x ,A y ,A z ), we obtain dn^ . . . 

—— = - Tn x + A x {p n + p 00 ) 

, at 

—^-=-Tn x + A x n + Sl y n z -Q z n y , (B6) + fl v {pu - Poo) - (B12) 

^=-rn, + A y n-^ + ^, (B7) ^ = - Tn, + A,( Pll + Poo) 

dn z -ttxipn - poo) + tl z n x , (B13) 

= - Tn z + A x n + f2 x % - Cl y n x , (B8) ^ 1 

, — 7— = - (r - A z )pu + -(A x n x + A y n y ) 

— = - 1 n + A x n x + A y n y + A z n z . (B9) ^ 

M +-(n x n y -n y n x ), (B14) 

In terms of the Bloch vector the qubit state population dpoo 1 
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[16] 



T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 [17 
(1998). 

E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lund- [18 
gren, and D. Preda, Science 292, 472 (2001). 
S. Suzuki and M. Okada, in Quantum Annealing and Re- [19 
lated Optimization Methods, edited by A. Das and B. K. 
Chakrabarti (Springer, 2005), vol. 679 of Lecture Notes [20 
in Physics, pp. 207 - 238. 

A. Das and B. K. Chakrabarti, Rev. Mod. Phys. 80, 1061 [21 
(2008). 

G. E. Santoro, R. Martonak, E. Tosatti, and R. Car, [22 
Science 295, 2427 (2002). [23 
M. Ohzeki and H. Nishimori, J. Comp. Theor. [24 
Nanoscience 8, 963 (2011). 

S. Tanaka and R. Tamura, in Lectures on Quantum Com- [25 
puting, Thermodynamics and Statistical Physics, edited 
by M. Nakahara and S. Tanaka (World Scientific, 2012), [26 
vol. 8 of Kinki University Series on Quantum Computing, 
pp. 3-62. " [27 

L. Stella, G. E. Santoro, and E. Tosatti, Phys. Rev. B 
72, 014303 (2005). [28 
S. Suzuki, H. Nishimori, and M. Suzuki, Phys. Rev. E 
75, 051112 (2007). [29 
M. H. S. Amin, Phys. Rev. Lett. 100, 130503 (2008). [30 
V. N. Smelyanskiy, U. v Toussaint, and D. A. Timucin, [31 
Dynamics of quantum adiabatic evolution algorithm for 
Number Partitioning, arXiv: quant-ph/0202155 [32 
T. Jorg, F. Krzakala, J. Kurchan, and A. C. Maggs, Phys. 
Rev. Lett. 101, 147204 (2008). [33 
A. P. Young, S. Knysh, and V. N. Smelyanskiy, Phys. 
Rev. Lett. 101, 170503 (2008). [34 
G. P. Berman and A. I. Nesterov, IJQI 7, 1469 (2009). 
A. I. Nesterov and G. P. Berman, Phys. Rev. A , (to [35 
appear in PRA) (2012). 

M. V. Berry, Proc. R. Soc. A 392, 45 (1984). [36 



M. V. Berry and M. Wilkinson, Proc. Roy. Soc. A 392, 
15 (1984). 

P. M. Morse and H. Feshbach, Methods of Theoretical 

Physics (McGraw-Hill, New York, 1953). 

A. P. Seyranian, O. N. Kirillov and A. A. Mailybaev, J. 

Phys. A: Math. Gen. 38, 1723 (2005). 

O. N. Kirillov, A. A. Mailybaev, and A. P. Seyranian, J. 

Phys. A 38, 5531 (2005). 

A. A. Mailybaev, O. N. Kirillov, and A. P. Seyranian, 

Doklady Math. 73, 129 (2006). 

S. Katsura, Phys. Rev. 127, 1508 (1962). 

J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005). 

A. C. M. Carollo and J. K. Pachos, Phys. Rev. Lett. 95, 

157203 (2005). 

A. I. Nesterov and S. G. Ovchinnikov, Phys. Rev. E 78, 
015202(R) (2008). 

E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407 
(1961). 

M. N. Lebedev, Special Functions & Their Applications 
(Dover, New York, 1972). 

L. Landau and E. M. Lifshitz, "Quantum Mechanics" 

(Pergamon, New York, 1958). 

C. Zener, Proc. R. Soc. A 137, 696 (1932). 

M. V. Berry, Proc. Roy. Soc. A 430, 405 (1990). 

A. Joye, G. Mileti, and C.-E. Pfister, Phys. Rev. A 44, 

4280 (1991). 

H. K. A. Joye and C. E. Pfister, Ann. Phys. 208, 299 
(1991). 

A. Kvitsinsky and S. Putterman, J. Math. Phys. 32, 1403 
(1991). 

R. Schilling, M. Vogelsberger, and D. A. Garanin, J. 
Phys. A 39, 13727 (2006). 

G. Dridi, S. Guerin, H. R. Jauslin, D. Viennot, and 
G. Jolicard, Phys. Rev. A 82, 022109 (2010). 
A. M. Dykhne, Sov. Phys.-JETP 14, 941 (1962). 



14 



[37] J. P. Davis and P. Pechukas, J. Chem. Phys. 64, 3129 

(1976) . 

[38] J.-T. Hwang and P. Pechukas, J. Chem. Phys. 67, 4640 

(1977) . 

[39] S. Suzuki and M. Okada, J. Phys. Soc. Japan 74, 1649 
(2005). 

[40] T. Caneva, R. Fazio, and G. E. Santoro, Phys. Rev. B 

76, 144427 (2007). 
[41] A. Erdelyi, W. Magnusand, and F. Oberhettinger, Higher 

Transcendental Functions, vol. I (McGraw-Hill, New 

York, 1953). 

[42] W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 

95, 105701 (pages 4) (2005). 
[43] J. Roland and N. J. Cerf, Phys. Rev. A 65, 042308 (2002). 



[44] M. Ohzeki, Phys. Rev. Lett. 105, 050401 (2010). 
[45] M. Ohzeki and H. Nishimori, Comp. Phys. Com. 182, 
257 (2011). 

[46] F. W. J. Olver, J. Res. Natl. Bur. Stand. B 63, 131 
(1959). 

[47] N. V. Vitanov and B. M. Garraway, Phys. Rev. A 53, 
4288 (1996). 

[48] Since Ylk>o cos ^ = 0> * ne term, Jcosfc, does not make 

contribution into the total Hamiltonian, and may be 

omitted in eo(t). 
[49] Approximation cos k « 1 instead of cos fe « 1 - k 2 /2 

leads to omission of corrections to the formula (|37p up to 

0(5/g). 



